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1 . 0  SUMMARY 

This  technical  report  presents  the  results  and  progress 
of  the  first  six  months  of  research  on  numerical  methods  for 
two-dimensional  process  modeling  under  Contract  MDA903-80-C-0498 , 
DARPA  Order  No.  3984.  A  significant  milestone  has  been  passed 
in  this  initial  period,  well  in  advance  of  our  anticipated 
schedule:  An  algorithm  of  sufficient  speed  to  be  used  for 
routine  VLSI  process  modeling  has  been  found.  It  is  orders 
of  magnitude  faster  than  methods  presently  being  used  in  the 
U.S.  and  abroad,  requiring  only  seconds  per  process  step  on 
the  CDC  Cyber  176  to  generate  full  two-dimensional  detail  on 
dopant  spread.  A  description  of  the  algorithm  will  be  found 
in  Section  3.2  below. 

1.1  TASK  OBJECTIVES 

The  overall  objective  of  this  program  is  to  develop  fast 
and  accurate  methods  for  computer  modeling  of  the  two- 
dimensional  spread  of  dopants  and  other  defects  during  VLSI 
circuit  fabrication.  Our  initial  goals  are  to  demonstrate 
these  methods  for  nonlinear  diffusion  of  a  single  dopant 
during  oxide  growth,  and  to  provide  the  resulting  computa¬ 
tional  algorithms  in  a  form  suitable  for  incorporation  into 
a  general  process  modeling  computer  code,  such  as  Stanford's 
program  SUPREM. 

Three  tasks  have  been  defined  for  the  first  year's  effort: 
Task  1  —  Analysis  of  Numerical  Methods 
Task  2  -  Formulation  of  Test  Cases 
Task  3  -  Algorithm  Development  and  Evaluation  . 
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Their  objectives  are,  respectively: 

Task  1  —  Select  promising  algorithms  from  other  disciplines, 
primarily  fluid  dynamics,  where  considerable 
progress  on  numerical  methods  for  multidimensional 
problems  has  recently  been  made. 

Task  2  —  a)  Determine  realistic  parameter  ranges  to  be 
used  as  test  cases  for  source  and  drain 
diffusion  into  a  short  MOSFET  channel  region. 


b)  Develop  and  analyze  two-dimensional  diffusion 
problems  for  which  dopant  profiles  can  be 
obtained  to  high  accuracy  by  existing 
techniques. 

Task  3  —  a)  Evaluate  each  selected  algorithm  for  speed 
and  accuracy. 

b)  Adapt  the  most  promising  algorithm (s)  to 
solve  the  combined  oxidation-nonlinear 
diffusion  problem  for  a  single  dopant 
species  in  two  dimensions. 

1.2  TECHNICAL  PROBLEM 

The  fabrication  of  VLSI  devices  requires  production  of 
features  of  submicron  size  and  separation.  Electrical  char¬ 
acteristics  such  as  threshold  and  punchthrough  voltages  will 
be  sensitive  to  dopant  spread  into  critical  areas  adjacent 
to  the  original  features.  Experimental  control  of  this 
spreading,  without  guidance  from  accurate  computer  modeling, 
will  be  costly,  tedious,  and  time-consuming.  However,  the 
use  of  standard  numerical  methods  to  achieve  an  adequate 
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modeling  capability  is  also  costly  and  time-consuming.  One 
should  therefore  seek  advanced  methods,  drawn  from  areas 
such  as  fluid  dynamics,  where  considerable  effort  and 
ingenuity  have  been  expended  in  recent  years  to  develop  fast 
and  accurate  solvers  for  the  characterization  of  multidimen¬ 
sional,  time-dependent  phenomena. 

1.3  GENERAL  METHODOLOGY 

Based  upon  our  own  ongoing  research  in  computational 
nonlinear  aerodynamics,  we  have  identified  several  promising 
approaches  to  the  development  of  a  fast  solver  for  two- 
dimensional  diffusion  problems.  After  a  preliminary  screen¬ 
ing,  a  few  of  these  have  been  selected  for  adaptation  to  the 
problem  of  dopant  spread  during  oxidation  or  annealing.  These 
algorithms  will  be  tested  for  speed  and  accuracy  on  the  prob¬ 
lem  of  nonlinear  dopant  diffusion  into  the  channel  region  of 
a  MOSFET,  as  well  as  on  simpler  problems  for  which  the  actual 
dopant  profiles  can  be  accurately  obtained  by  other  means. 

1 . 4  TECHNICAL  RESULTS 

Two  algorithms  have  gone  through  preliminary  screening: 
high-order  finite  element  methods,  and  a  multidimensional 
version  of  the  method  of  lines,  both  utilizing  an  optimized 
stiff  integrator  for  the  time  integration.  The  finite  element 
methods  have  proved  disappointing,  but  the  method  of  lines 
has  provided  an  unexpectedly  large  gain  in  speed.  Two- 
dimensional  post-annealing  profiles  (without  oxidation)  for 
a  15-minute  drive-in  cycle  on  boron-implanted  silicon  at 
1100°C,  with  a  peak  initial  boron  concentration  of  1020  cm-3, 
were  obtained  in  ten  seconds  of  CPU  time  on  the  Cyber  176. 
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These  results  can  be  directly  compared  with  the  recently 
published  work  of  Warner  and  Wilson  (Bell  System  Tech.  J.  59 , 
1  (1980)),  where  similar  nonlinear  diffusion  problems  with 
the  same  number  of  unknowns  (a  21 x  41  grid) ,  solved  by 
second-order  finite  element  methods,  took  over  seven  minutes 
on  the  Cray-1  computer,  a  machine  which  is  five  to  fifteen 
times  faster  than  the  Cyber  176.  Comparing  equivalent  CPU 
times  on  the  Cyber  176,  one  finds  our  method  of  lines 
algorithm  to  be  at  least  two  orders  of  magnitude  faster 
than  the  algorithm  used  by  Warner  and  Wilson.  A  more 
detailed  discussion  of  these  results  can  be  found  in 
Sections  3.2  and  3.3. 

1.5  IMPORTANT  FINDINGS  AND  CONCLUSIONS 

The  speed  of  the  method  of  lines  algorithm  is  already 
sufficient  to  render  its  use  for  routine  process  design 
practical  for  a  process  engineer  with  access  to  the  latest 
IBM  or  CDC  computers.  Since  these  machines  can  currently 
be  rented  economically  on  a  time-sharing  basis,  it  appears 
that  we  now  have  in  hand  the  core  of  a  new  two-dimensional 
process  simulator.  Our  immediate  priorities  should  therefore 
be  to  complete  the  testing  and  characterization  of  this 
algorithm,  and  to  extend  its  use  to  diffusion  during 
oxidizing  process  cycles,  where  the  oxide-silicon  boundary 
moves  nonuni formly. 

1.6  SPECIAL  COMMENTS 

To  facilitate  the  transfer  of  this  algorithm  into  a 
complete,  two-dimensional  process  simulator,  an  early  meeting 
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with  researchers  at  Stanford  responsible  for  simulator 
development  has  been  tentatively  scheduled  for  March,  1981. 

By  this  date,  good  measures  of  the  overall  performance  of 
the  method  of  lines  algorithm  will  have  been  obtained. 

1.7  IMPLICATIONS  FOR  FURTHER  RESEARCH 

The  work  of  testing  potentially  fast  algorithms  for 
multidimensional  VLSI  problems  has  barely  begun.  There 
remain  at  least  two  strong  motivations  to  continue  research 
on  new  methods:  to  find  methods  which  can  be  applied 
economically  to  derive  VLSI  device  electrical  characteristics 
from  two-dimensional  dopant  profiles,  and  to  extend  the 
treatment  of  device  structure  formation  to  the  level  of 
the  actual  defect  processes,  involving  multispecies  migration 
and  chemical  reaction.  This  will  ultimately  provide  a  better 
bridge  from  the  physics  of  solid-state  defect  dynamics  to 
the  engineering  observables  than  the  current,  semi-empirical 
models.  Such  a  tool  can  be  used  to  recalibrate  these  models 

for  new  situations,  and  to  determine  at  what  scale  new  j 

i 

effects  will  appear  and  become  important  to  device  performance. 

1.8  ORGANIZATION  OF  THE  REPORT 

The  basic  formulation  of  the  equations  describing  one¬ 
dimensional  nonlinear  diffusion,  with  segregation  at  a  moving 
oxide-silicon  boundary,  is  given  in  Section  2.1.  This  is 
followed  by  a  presentation  of  the  two  numerical  procedures 
examined  to  date  {Section  2.2)  and  a  summary  of  results  for 
one-dimensional  test  cases  (Section  2.3). 
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In  Section  3,  the  generalization  of  these  problems  to 
two  dimensions  is  considered.  For  the  reasons  stated  in 
Section  3.2,  two-dimensional  numerical  tests  have  only  been 
carried  out  on  the  method  of  lines  algorithm.  Alternative 
approaches  to  the  treatment  of  the  nonuni formly  moving  oxide- 
silicon  boundary  are  described  in  Section  3.1.  Results  for 
a  nonmoving,  planar  boundary  are  presented  in  Section  3.3. 
Plans  for  future  work  are  outlined  in  Section  4. 
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2.0  ONE-DIMENSIONAL  REDISTRIBUTION  PROBLEMS 

Initial  computations  were  carried  out  on  one-dimensional 
redistribution  problems  to  keep  computer  costs  low,  to  learn 
the  limitations  of  the  various  numerical  procedures,  and  to 
determine  the  number  of  grid  points  required  to  obtain 
satisfactory  resolution. 

2.1  MATHEMATICAL  FORMULATION 

The  redistribution  of  impurities  in  silicon-on-sapphire 
(SOS) -type  structures  during  thermal  oxidation  was  treated 
by  Maldonado  and  Murphy^  as  a  nonlinear  diffusion  equation 
with  a  moving  boundary.  This  equation  can  be  transformed  to 
one  having  a  fixed  boundary  by  introducing  the  change  of 
variables 

^  =  (x-mU(t)  )LI/L(t)  (la) 

x  =.t  (lb) 

where  m  is  the  ratio  of  the  thickness  of  silicon  consumed, 
given  by  L(0)  -L(t),  to  the  oxide  thickness,  U(t).  L^  repre¬ 
sents  a  characteristic  length  at  t=0,  which  is  given  by 
LI  =  L(0)  -mU(0)  and  denotes  the  right  boundary.  In  these 
new  variables  the  governing  equation  for  the  impurity 
concentration,  N(£,x),  is  given  by 
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with  boundary  conditions 


3N  (0  ,  t)  =  L  (t) 
3x  Li 

3N(Lt,t) 

1  _ 


(k  -  m) 


U  (t  )  N  ( 0 ,  t ) 
D  (N  (0 , t)  ) 


(2b) 

(2c) 


and  initial  condition 


N(?,0)  =  N°  (£)  .  (2d) 

k  is  the  segregation  coefficient,  and  N°  (£)  is  usually  selected  as 
a  Gaussian,  modeling  the  type  of  doping  profile  which  results 
from  ion  implantation.  The  nonlinear  diffusion  function  is1 


D  (N ) 


N  + 


N2  +  2n? 
_ in 

(N2  +  4n?n) ** 


(2e) 


where  n^n  is  the  intrinsic  carrier  density  at  the  given 
process  temperature  and  D°  is  the  intrinsic  diffusivity. 

This  nonmoving  boundary  value  problem  (2)  is  much  easier 
to  solve  numerically  than  the  original  one  proposed  by 
Maldonado  and  Murphy.^-  In  fact,  the  numerical  results  to  be 
presented  will  indicate  that  good  accuracy  is  achievable 
using  relatively  coarse  grids  when  discretizing  equation  (2) . 
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2.2  NUMERICAL  PROCEDURES 


2.2.1  B-splines 


The  Galerkin  method  using  piecewise  polynomial  B-splines 
2  . 

(see  deBoor  )  consists  of  approximating  the  concentration  by 


N  (?,t)  =  £  a. (t)B. (?) 
s  i=l  1  1 


where  the  B^(?)  are  known  polynomials  of  order  k  which  vanish 
outside  of  specified  intervals.  If  N  is  replaced  by  Ng  in 
equation  (2)  and  the  resulting  expression  is  multiplied  by 
Bq(?)  and  integrated  over  [°»Lj  >  the  Galerkin  formulation 
takes  the  form 


fLi  3Ns(?,t) 

J  _ 


LI  l2  f4 

-  nryj  J0  D 


Bq(?)ds  = 


3N  dB 

<V  TT  d#  d5 


I*  .  rL  I  3N  L 

+  tjtt  J0  (Li-e)  -5T  Vc  •  nrr  o'-'”>u<,)Ns«>,T>Bq(o> 

(4a) 

fLl  rLl 

J  N  ( ? ,  0 ) B  (? )  d?  =  I  N° (?) B  (?) d?  (4b) 


q  =  1, 2, • • • ,n 
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where  one  integration  by  parts  was  applied.  Equation  (4)  may 
be  written  in  matrix  form  as  the  following  system  of  nonlinear 
initial  value  ordinary  differential  equations: 


F (a, t) 


(5) 


where  a  =  (a^(r),  o^ff) 


A  =  | fo  \<UB} 


(5)dC|, 


T 

*  •  • ,  otn  (t)  )  , 

and  the  initial  value  a(0) 


is  obtained 


from  solving  equation  (4b) . 

The  following  numerical  characteristics  are  well  known 
for  such  problems  as  (4)  and  (5) : 


(1)  The  matrix  A  has  half  bandwidth  k-1. 

(2)  All  integrals  in  equation  (4)  may  be  evaluated 

1. 1. 

accurately  using  k  order  Gaussian  quadratures 
between  the  knots  of  the  B-splines  if  the  separation 
is  sufficiently  small. 

(3)  Equation  (5)  is  stiff  (see  Gear1 2 3 4)  and,  consequently, 
stiffly  stable  integration  procedures  must  be  used. 
This  typically  requires  the  evaluation  of  the 
Jacobian,  8F/3a,  which  also  has  half  bandwidth 
equal  to  k-1. 

(4)  k  order  B-splines  lead  to  numerical  approxima- 

th 

tions,  Ng,  which  are  k  order  accurate  in  the 
spatial  discretization  parameter. 


The  time  integration  of  (5)  is  carried  out  using  a 
robust  integrator  with  automatic  variable  step-and-order 
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changes  from  first  to  fifth  order.  Thus,  the  user  need  only 
specify  an  error  tolerance,  and  the  integrator  automatically 
guarantees  that  the  time-integration  steps  satisfy  this 
criterion.  Spatial  errors  are  reduced  by  increasing  the 
number  n  of  B-splines  or  increasing  their  order,  k.  The 
motivation  for  considering  such  a  method  is  that  by  employing 
high  order  B-splines  (large  values  of  k)  we  might  be  able  to 
reduce  the  number  n  required  to  provide  acceptable  accuracy 
and  as  a  by-product  reduce  computer  time.  This  must  be 
balanced  against  any  increase  in  computation  caused  by  the 
increase  in  bandwidth  of  the  matrices  A  and  9F/3ot. 

2.2.2  Method  of  Lines 

We  will  describe  this  method  by  using  a  uniform  mesh 
A£  =  Lj/n  on  the  interval  although  a  nonuniform  mesh 

has  been  applied,  occasionally.  The  partial  differential 
equation  (2a)  is  discretized  in  the  spatial  variable  £  using 
centered  differences;  i.e.,  for  q  =  l,2,***,n-l. 


+ 


mU  (t) 
L(t) 


(6a) 


where  Cq  =  qA£,  Nq  =  N(5q,T>,  Kq±h  =  |(Nq±1  +  Nq),  and 
D  +.  =  D (N  +J) .  The  boundary  conditions  are  discretized 
using  one-sided  differences  in  (2a) ,  and  combining  this 
with  (2b)  and  (2c)  yields 
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inverted  easily.  This  integrator  is  also  variable-order  and 
variable-step-size  and  employs  the  implicit  backward  differ¬ 
entiation  formulas,  which  are  now  quite  common  for  treating 
stiff  ordinary  differential  equations.  Excellent  error 
control  is  the  strong  point  for  this  technique. 

2 . 3  NUMERICAL  RESULTS 

I 

i 

2.3.1  B-splines 

For  our  test  case  we  use  an  example  described  by  Prince 
and  Schwettmann. 5  This  involves  the  redistribution  of  a  high 
dose  80  keV  boron  implant  in  (111)  silicon.  The  redistribution 
is  done  by  bringing  the  silicon  into  contact  with  an  oxidizing 
(steam)  ambient  at  1100°C.  This  is  a  highly  nonlinear  case 
and  contains  a  highly  peaked  Gaussian  as  an  initial  condition, 
which  is  difficult  to  approximate  by  an  equation  of  the  form 
(3)  when  n  is  relatively  small.  To  alleviate  this  problem  we 
employ  a  logarithmic  transformation  of  the  form 

N(5,t)  =  eu(C'T)  . 

In  Table  1  below,  a  list  is  given  of  some  solution  values  at 
the  Si02~Si  interface  at  time  t=0.75  hour  for  various  order 
k  and  number  of  uniform  subintervals  it  on  a  silicon  slab  of 
length  2  ym.  Here  n  =  Jt  +  k  -  1.  CPU  denotes  the  running 
time  on  the  IBM  3033  in  minutes  and  3.169(18)  signifies  the 
number  3.169  *  1018 . 

Observe  that  the  CPU  time  is  approximately  linear  in  £ 
but  quadratic  in  k.  Also,  since  the  solution  ranges  over 
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Table  1.  Solution  Using  B-splines 


k 

SL 

CPU 

u 

N 

2 

50 

0.089 

42.60 

3.169 (18) 

3 

20 

0.080 

42.70 

3.503(18) 

3 

30 

0.117 

42.58 

3.106  (18) 

3 

50 

0.191 

42.57 

3.076  (18) 

4 

10 

0.061 

43.17 

5.604  (18) 

4 

25 

0.222 

42.61 

3.201(18) 

4 

50 

0.400 

42.57 

3.076  (18) 

5 

50 

0.682 

42.57 

3.076  (18) 

6 

25 

0.578 

42.58 

3.106  (18) 

6 

50 

1.081 

42.57 

3.076(18) 

many  orders  of  magnitude,  a  certain  minimum  number  of  subintervals 
is  required  to  accurately  monitor  the  redistribution.  In  this 
example  the  third  order  method  (k=3)  employing  30  subintervals 
is  about  optimal  when  both  accuracy  and  computer  time  are 
considered.  This  means  that  in  two  dimensions  a  grid  of 
30  x  60  points  using  a  third  order  method  should  suffice  when 
the  lateral  dimension  is  twice  the  depth  (a  typical  applica¬ 
tion  in  MOSFET  process  design) . 

2.3.2  Method  of  Lines 

The  logarithm  transformation  needed  for  high-order 
B-splines  was  not  required  here,  since  the  initial  data  (7d) 
is  used  directly  as  initial  conditions  for  the  ordinary 
differential  equations.  In  Table  2  we  again  list  some 
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Table  2.  Solution  By  Method  of  Lines 


n+1 

CPU 

N 

205 

0.347 

6.039  (18) 

45 

0.107 

6.040 (18) 

21 

0.068 

6.152(18) 

interface  values  at  time  x  =  0.75  hour  for  various  numbers 
of  grid  points  n+1.  However,  this  time  we  use  the  SUPREM 
default  values  for  the  parameters  in  the  code,  which  explains 
the  differences  between  the  surface  concentration  values  in 
the  two  tables. 

SCJPREM  gave  a  value  of  N  =  5.49(18)  using  0.076  minute 
of  CPU  time.  Observe  that  acceptable  accuracy  can  be  achieved 
with  only  21  grid  points,  which  translates  to  a  grid  of 
21 x  41  in  two  dimensions. 
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3.0  TWO-DIMENSIONAL  REDISTRIBUTION  PROBLEMS 

As  we  have  learned  by  studying  one-dimensional  problems, 
the  moving  boundary  problem  can  be  greatly  simplified  by 
mapping  it  into  a  fixed  region.  In  two  dimensions,  two 
different  types  of  mappings  into  a  rectangle  are  considered: 
(1)  translation-stretching,  and  (2)  conformal  mapping. 

3.1  FORMULATION 

3.1.1  Translation-Stretching 

In  Figure  1  below  is  a  model  for  the  field  oxidation 
step  in  the  fabrication  of  an  MOS  device  in  physical 
coordinates. 


SC80  11378 


Figure  1.  Definition  of  Symbols  for  Three-Region  Model 
of  Field  Oxidation 
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In  two  dimensions  the  moving  boundary  is  characterized 
by  the  function  mU(y,t).  The  above  region  is  mapped  into  a 
rectangle  by  the  transformation 


x  -  mU  (y ,  t) 


L  (y,  t) 


n  =  y 


T  =  t 


where  L(y,t)  =  L0-mU{y,t),  Lj  =  LQ-mU(y,0)  =  Lg-mU^,  and 
U(y,0)  =  UQI  =  Uqjj  =  uoiii  =  constant.  Under  this  transforma¬ 
tion  the  standard  two-dimensional  nonlinear  diffusion  equation 
can  be  written  as 


Li  +  m¥(LT  -  5)21  a  a  r  aw 

“ - ^ -  &  [D(N)  IfJ  +  ^lD(N)  I? 


2mU  a 

V1  «•!-«>  &  D(N) 


m(L  -  C) 

+  - - -  U  L-(U  L  +  2mU2  )  D  (N) 

L2  l  t  hn  n 


9N 

35 


for  0  <  £  <  Lj  ,  O^ti^b^ 


with  boundary  conditions 
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9N(5,0,t)  _  „ 

an  "  0 


0  <  5  <  L. 


(9b) 


3N(5,b3,x) 


an 


=  o 


0  <  5  <  L, 


(9c) 


3N(LI#n,T) 

95 


=  0 


0  <  n  ^  b. 


(9d) 


LjD (N)  ,N 

L  3{ 


fmtl 


5=0 


IT  LID(N)  H 


5=o 


U  D  (N) 

n 


5=o 


=  (U2  +  l)(k-m)UxN(0,rifT)  ;  0  <  n  <  b3  . 


(9e) 


The  nonlinear  diffusion  coefficient  is  selected  as  that 

7 

used  by  Warner  and  Wilson: 


D(y) 


jl  +  19  [y  +  \/y2  +  l]j 

1  +  -  ^ 

(  ”  ) 

[  v/Y2  +  lJ 

( 9  f ) 


where  y  =  N/2n^n* 

As  a  test  case  to  characterize  the  moving  boundary 
function  U(n»t),  define 


Ui(t)  =  -  (a/2)  +  [a2/4  +  St]*5 

(10a) 

UrII  (  t )  =  3ur(t) 

(10b) 
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where  a  and  0  are  oxide  growth  parameters  in  the  standard 
Deal-Grove  model.  To  interpolate  smoothly  between  these 
limits,  take 

u(n,T)  =  ux(t)  ;  0  <  n  <  bx  (10c) 


<uiii"ui>  <um 

U(n'T)  *  UI  +  . 


-U  )  [(n-b-Jir 

- — —  sin  2  -r- - c — 

11  t>2  —  fc>i 


b1  <  n  <  b2  (lOd) 


u(n,T)  =  uII3;(t)  ;  b2  <  n  <  b3 


(lOe) 


Although  many  other  selections  are  possible,  this  one  has 
the  advantage  that  U  is  a  continuous  function  with  continuous 
first  derivative,  and  that  equation  (9a)  is  greatly  simplified 
in  the  regions  0  <  n  <  b^  and  b2  <  n  <  b3  where  =  0. 


3.1.2  Conformal  Mat 


The  introduction  of  mixed  derivative  terms  02/a£9n) 


into  the  diffusion  equation  and  the  presence  of  both  and 
n-derivatives  in  the  boundary  condition  at  the  oxide-silicon 
interface  are  consequences  of  the  lack  of  orthogonality 
between  lines  of  constant  5  and  lines  of  constant  ti  in  the 


translation-stretching  transformation.  The  impact  of  these 
complications  on  the  efficiency  of  the  numerical  procedure 
is  not  small;  taken  together,  they  have  been  found  to  nearly 
triple  the  CPU  time  required  to  obtain  a  solution. 
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A  simple  approach  to  overcome  these  difficulties  is  to 
employ  a  conformal  map  whose  level  lines  closely  approximate 
the  desired  shape  of  the  silicon  region.  The  construction 
of  appropriate  conformal  maps  for  the  area  beneath  a  bird's 
beak  oxide  layer  was  studied  during  FY  1980  under  Rockwell 
IR&D  funding;  the  application  of  one  such  mapping  to  this 
problem  is  summarized  below. 

One  can  find  in  standard  texts  on  complex  functions^  the 
transformation 


z  =  cosh-1  (2-Tr--r 


1  \ 


(k  +  1)  5  -  2k 


(k-  l)c 


Im  c  >  0  ,  (11) 


which  for  k  >  1  maps  the  upper  half  £-plane  onto  a  strip  with 
a  step  in  it: 


tt  (k  -  l)/k  <  y  <  it  for  x  <  0  ,  and  0  <  y  <  tt  for  x  >  0  . 

Since  the  exponential  transformation  £  =  ew  maps  a  uniform 
strip  of  width  tt  (0  <Imw<  tt)  onto  the  upper  half  £ -plane, 
the  indicated  substitution  of  ew  for  z;  in  (11)  yields  a 
conformal  map  from  a  uniform  strip  to  a  stepped  strip.  The 
exact  correspondence  is  illustrated  in  Figure  2.  A  level 
line  Imw=  >  0  provides  a  natural  model  for  the  oxide- 
silicon  boundary,  the  region  above  £q  corresponding  to  the 
silicon  layer. 

The  basic  advantage  of  employing  conformal  maps  is  that 
the  structure  of  neither  the  differential  operator  nor  its 
associated  boundary  conditions  is  changed  thereby.  The 
diffusion  equation 
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Figure  2.  Conformal  Map  from  Strip  to  Stepped  Strip: 
w=n  +  i£»  z  =  x+iy 


(12) 


is  transformed  into 


dw  2jiL/n  +  _L/n  =  iN  +  iN3£  3N9n 

dz  |9£\  9U  3n\  9n/  3t  9£  9t  9n  9t  ' 


(13) 


where  the  time  dependence  of  £  and  n  arise  from  motion  of 
the  oxide-silicon  boundary.  Using  (11) ,  the  functions 
9£/9t,  9d/9t,  and  dw/dz  are  expressible  directly  in  terms  of 
n,  and  t,  so  that  computations  of  dopant  diffusion  during 
oxidation  can  be  carried  out  on  a  constant  (£,ri)  grid  without 
reference  to  physical  coordinates.  A  single  transformation 
of  the  dopant  profiles  to  physical  coordinates  can  be 
performed  after  the  process  cycle  is  complete. 
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The  boundary  conditions  at  the  top  and  bottom  of  the 
silicon  layer  are  straightforward  to  derive,  and  involve  only 
the  normal  derivative  3N/3£.  However,  the  conditions  at  the 
lateral  edges  of  the  computational  domain  require  some  care. 

If  lines  | n |  =  constant  are  chosen  to  bound  this  domain, 
then  the  physical  boundaries  in  x  and  y  will  be  somewhat 
curved,  and  will  change  position  as  the  oxide-silicon  boundary 
moves.  Fortunately,  the  effect  of  conditions  at  these 
boundaries  on  the  dopant  profiles  in  the  region  of  interest 
can  be  made  small  by  expanding  the  computational  domain 
laterally.  More  accurate  procedures  are  available,  but  are 
not  needed  in  the  present  context. 

The  basic  question  to  be  answered  with  regard  to  the 
use  of  conformal  maps  for  this  problem  is  how  large  a 
computational  burden  is  imposed  by  the  necessity  to  evaluate 
the  functions  3£/3x,  3h/9t,  and  dw/dz  at  each  time  step. 

This  must  be  compared  against  the  approximate  tripling  of 
CPU  time  experienced  with  the  translation-stretching 
transformation  for  the  nonuni formly  moving  boundary  problem. 

3.2  NUMERICAL  PROCEDURES  IN  TWO  SPATIAL  DIMENSIONS 

3.2.1  B-splines 

Although  the  method  of  B-splines  may  be  extended  to 
two  dimensions,  a  number  of  disadvantages  make  this  method 
noncompetitive  with  the  method  of  lines.  The  first  difficulty 
is  large  storage  requirements.  For  example,  using  our 
estimate  of  a  grid  of  30  *  60  with  a  third  order  B-spline  in 
each  of  the  two  directions,  the  matrices  A  (equation  (5)) 
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and  3F/3a  each  have  half  bandwidth  w  =  2  +  2(30  +  2)  =  66. 

Because  of  pivoting  requirements,  the  total  storage  for  both 
matrices  A  and  3F/3a  is  s  =  2  *  32  *  62  *  (3w+l)  =  789,632  words. 

If  the  smaller  20  *  40  grid  with  second  order  B-splines  were 
adequate,  then  w  =  21  and  s  =  2*21*41*  (3w+l)  =  110,208  words. 
In  either  case  the  storage  requirements  are  excessive. 

Secondly,  to  invert  these  banded  matrices  with  relatively 
large  bandwidth  many  times  would  require  an  excessive  amount 

7 

of  computer  time.  In  fact,  Warner  and  Wilson  have  solved 
the  easier  nonmoving  boundary  value  problem  using  second 
order  B-splines  on  a  21  *  41  grid  and  their  worst  cases 
required  7.25  minutes  on  a  Cray-1  computer  to  integrate  to 
t  =  0.25.  This  translates  to  at  least  36.25  minutes  on  a 
CDC  176  or  an  IBM  3033*.  Finally,  the  cross  derivative  term 
in  equation  (9a)  may  lead  to  some  difficulty  because  the 
existence  of  Galerkin  solutions  to  such  problems  have  not 
completely  been  determined.  Of  course,  the  conformal  mapping 
formulation  doesn't  suffer  from  this  disadvantage,  but  large 
storage  requirements  still  limit  the  B-spline  applicability. 


*Time  comparisons  between  conventional  mainframe  computers 
and  the  Cray-1  are  complicated  by  the  pipeline  architecture 
of  the  Cray,  which  can  be  used  to  advantage  on  large  vector 
problems.  The  quoted  factor  of  five  above  assumes  that  no 
modification  of  the  original  FORTRAN  code  to  utilize  this 
feature  was  attempted  by  Warner  and  Wilson. 
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3.2.2  Method  of  Lines 

Consider  a  uniform  mesh  in  the  C  and  n  direction  with 
meshwidth  AC  and  An,  respectively.  For  interior  points  the 
spatial  derivatives  in  equation  (9)  are  discretized  using 
the  following  centered  difference  approximations: 


9N . 

_ ii 

3C 


i 


N 


-  N. 


i+l.j  i-1 , j 
—  iKtr — 


(14a) 


3 

3C 


D(N) 


i+*5,  j 


(N.  ,  ,  .  -  N.  .) 

(AC)2 


j 


(N.  .  -  N.  .  .) 

*1  3--1/3 

(AC)  2 


(14b) 


Di  i  j+3s 


(N. 

1,1+1 


-  N.  .) 
H 


(An) 


(N.  .  -  N.  •  .j  ) 

D.  .  _±J - LUZi. 

1,1-H 


(An) 


(14c) 


3n1 

U] 


Di,  j+l(Ni+l,  j+1  ~  Ni-1,  j  +  l}  ~  Pi,  j-l(Ni+l,  j-1  ~Ni-l,j-l} 


4ACAn 


(14d) 


where  =  N(Ci,hj»x),  =  D  (N^  j )  *  etc. 

Including  the  boundary  conditions,  especially  equation  (9e) , 
leads  to  very  complicated  difference  expressions.  We  will 
merely  illustrate  the  basic  idea  using  one  of  the  second 
partial  derivative  terms.  Write  the  first  order  approximation 
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9 

n 


-  V£ 


AC 


(14e) 


where  the  term  9N1  ./3C  is  determined  from  equation  (9e) . 

In  the  regions  where  =  0,  equation  (9e)  is  much  simpler 
and  a  discrete  version  can  be  written  as 


9N  _  Ni+l,j  "  Ni-l,j 

ac  ~  2ac 


(k  -  m)U  N 


(14f) 


The  last  expression  may  be  solved  for  N.  .  .  and  the  result 

9  J 

is  substituted  into  equation  (14b)  to  produce  a  second  order 
approximation  in  the  regions  0  <  n  <  b^  and  b2  <  n  <  b3< 
Similar  techniques  are  used  for  the  other  boundary  conditions 
and  partial  derivatives. 

As  in  the  one-dimensional  case,  the  spatial  variable 
differencing  leads  to  a  semidiscrete  system  of  nonlinear 
ordinary  differential  equations.  The  equations  corresponding 
to  the  interior  mesh  points  have  the  form 


dN. 


J-'J  =  f.  .(N.  N.  ,  N .  ,  .  N.  »  .  .  ,  N.  .  , 
dx  if]  i,j  i-l,j  l+l,  j  1,3+1  i,j-l 


Ni+1, j+1'  Ni-1, j+1'  Ni+1, j-1'  Ni-l,j-lf  T)  *  (15) 


Similar  equations  may  be  written  for  the  boundary  points. 
Unfortunately,  equation  (15)  is  stiff,  and  consequently, 
the  Jacobian,  9f/9N,  is  needed  to  converge  the  corrector 
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equation  in  the  linear  multistep  method  used  to  solve  it. 

The  coefficient  matrix  in  this  Newton-like  method  is 

P  =  I  -  ATBQaf/8N  (16) 

where  Bq  is  a  scalar  associated  with  the  order  of  the  corrector 
equation.  Most  of  the  computer  time  is  used  in  solving  the 
linear  system 

Px  =  b  .  (17) 

Therefore,  it  is  critical  that  efficient  methods  be  developed 
for  solving  equation  (17) .  Fortunately,  the  matrix  P  need 
not  be  exact  but  may  be  approximated  by  a  function  which 
contains  the  "basic  nature"  of  the  problem  being  solved. 

The  consequence  of  using  an  approximate  Jacobian  is  that 
the  iterative  procedure  may  take  longer  to  converge,  but 
the  answer  will  still  be  correct.  This  is  analogous  to  using 
the  secant  method  in  place  of  Newton's  method  when  solving 
nonlinear  equations. 

Two  simplifications  are  made  in  approximating  the  Jacobian. 
The  first  is  to  ignore  the  cross-derivative  term  equation  (14d) 
when  evaluating  3f/9N  but  not  when  evaluating  f  (equation  (15)) . 
The  rationale  in  doing  this  is  that  the  cross-derivative  term 
is  only  present  in  a  very  small  region  of  the  computational 
domain  (b2  <  n  <  b^)  where  0,  so  that  its  absence  should 

not  greatly  affect  the  convergence  rate.  On  the  other  hand, 
if  this  term  were  present,  the  solution  of  equation  (17)  would 
have  to  be  done  by  the  relatively  slow  banded  matrix  techniques 
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instead  of  the  much  faster  iterative  methods.  The  second 
simplification  assumes  that  on  the  boundary  £  =  0  and 
<  n  <  b2  we  may  write  9D/3N  *  0.  There  is  no  real 
justification  for  making  this  boundary  approximation  other 
than  it  simplifies  the  evaluation  of  the  Jacobian.  The 
price  paid  for  such  an  approximation  is  slower  convergence. 

Applying  the  above  approximations  results  in  a  Jacobian 
with  only  5  nonzero  diagonals.  Successive  overrelaxation 
(SOR)  methods  may  be  employed  to  solve  equation  (17) .  The 
storage  problem  is  minimal,  and  as  P  changes  from  time  step 
to  time  step  very  accurate  initial  solutions  provide  fast 
convergence.  In  contrast,  banded  matrix  solvers  must  start 
from  scratch  each  time  step  and  do  not  use  any  previous 
information. 

The  numerical  solution  of  (15)  is  performed  by  a 
variable-order  and  variable-step-size  stiff  integrator 
similar  to  the  one-dimensional  integrator  with  the  critical 
equation  (17)  now  being  solved  by  SOR  rather  than  banded 
matrix  techniques.  Our  test  cases  show  that  this  procedure 
is  superior  to  the  B-spline  approach. 


3.3  NUMERICAL  RESULTS  IN  TWO  DIMENSIONS 

The  initial  condition  (ion  implant)  for  the  following 
test  cases  is 


N (x,y,0)  = 


Nd  x  lO4 
— — c —  exp 
(2tt) 


Nd  x  lO* 


<*- V* 


+  Nb  ;  0  <  y  <  a 
0  <  x  <  L, 


(2ir)  a. 


^ —  exp  pJ  (y  -  a)  2j  exp 


<*  -  v 


2°p  J 


+  Nv 


b  ' 

a  <  y  <  b,,  0  <  x  <  Ln 
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where  a  =  0.075  ym,  p  =  \l~2/o  ,  N,  =  5x  1014  cm-3,  a  =  b,/2, 

P  P  ^ 

N,  =  2x10 15  cm-2,  and  R  =  0.285  ym. 
d  p 

In  order  to  determine  the  lower  bounds  on  computer  time, 
we  first  consider  a  nonmoving  boundary  value  problem  studied 

7 

by  Warner  and  Wilson.  As  noted  earlier,  their  most  difficult 
case  required  about  36.25  minutes  on  a  CDC  176  to  integrate 
from  t  =  0  to  t=0.25.  Here  Lq  =  2.5,  b^  =  5 ,  UT  =  U  =  0, 

D°  =  0.0589  ym2/hr,  and  n^ =  9 . 25 x  io 18  cm-3 .  It  should  be 
noted  that  the  Jacobian  matrix  for  this  problem  is  exact 
since  no  cross  derivatives  are  present  in  equation  (9)  and 
all  boundary  conditions  are  homogeneous.  Therefore,  our 
computer  times  should  be  a  minimum  for  this  case.  The 
numerical  results  are  summarized  in  Table  3  for  two  computa¬ 
tional  grids.  Here  NSTEP  denotes  the  total  number  of  time 
steps  to  integrate  equation  (15) .  NFE  is  the  number  of 
functional  evaluations  required  of  the  integrator,  i.e., 
number  of  times  f  in  equation  (15)  is  evaluated.  NJE  is 
the  number  of  times  the  Jacobian,  9f/9N,  is  evaluated. 

Nil  is  the  total  number  of  internal  iterations  needed  to 
solve  equation  (17)  by  SOR  for  all  time  steps  up  to  t. 

NQ  denotes  the  present  order  of  the  integrator  (1  <  NQ  <  5) . 

At  is  the  present  step  size,  and  CPU  denotes  the  total  CPU 
time  in  seconds  to  solve  the  problem  on  a  CDC  176.  These 
times  are  remarkably  fast  compared  to  those  reported  by 
other  researchers.  Note  that  the  smallest  time  steps  occur 
between  0  <  x  <  0.05,  very  few  Jacobian  evaluations  are 
needed,  and  the  total  number  of  internal  iterations  (Nil) 
is  relatively  small. 

The  initial,  as-implemented  dopant  profile  and  the 
result  of  the  quarter-hour  drive-in  step  are  shown  in 
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Table  3  (a) .  11 x  21  Grid 


T 

NSTEP 

NFE 

NJE 

Nil 

NQ 

At 

CPU  (sec) 

0.05 

43 

7 

77 

4 

0.0055 

0.10 

55 

8 

100 

4 

0.0085 

0.15 

44 

63 

9 

115 

4 

0.0125 

0.20 

48 

67 

10 

123 

4 

0.0192 

0.25 

50 

71 

- 

10 

131 

4 

0.0192 

1.807 

Table 

3(b)  . 

21  x 

J1  Grid 

T 

NSTEP 

NFE 

NJE 

Nil 

NQ 

At 

CPU  (sec) 

0.05 

39 

59 

8 

4 

0.0040 

0.10 

48 

70 

10 

4 

0.0081 

0.15 

54 

77 

10 

If! 

5 

0.0113 

0.20 

58 

81 

10 

153 

5 

0.0113 

0.25 

62 

86 

11 

165 

5 

0.0200 

8.217 

Figures  3  and  4,  respectively.  The  lateral  spread  of  the  boron 
is  quite  evident,  amounting  to  a  few  tenths  of  a  micron. 
Accurate  information  of  this  kind  is  essential  to  prediction 
and  control  of  device  electrical  characteristics  when  the 
overall  channel  length  (separation  of  source  and  drain)  is 
of  the  order  of  one  micron. 

For  our  next  degree  of  difficulty  we  include  a  non- 
uniformly  moving  boundary  with  a  =  0.0914  and  3  =  0.576 
(equation  (10a)),  but  impose  Dirichlet  boundary  conditions 
at  5  =  0.  This  will  allow  us  to  study  the  effect  of  leaving 
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out  the  cross-derivative  term  in  the  evaluation  of  the 
Jacobian.  However,  the  Jacobian  of  the  boundary  conditions 
is  still  exact.  The  cross-derivative  approximation  equation 
(14d)  requires  a  finer  grid,  so  we  set  LQ=1.25  and  b3=2.5 
and  use  a  21  *  41  grid.  Computatioral  results  are  given  in 
Table  4.  In  comparing  Table  3(b)  with  Table  4,  we  note  that 
the  big  difference  is  the  number  of  internal  iterations, 

Nil,  which  is  larger  when  an  approximate  Jacobian  is  employed. 
The  integration  was  carried  out  to  a  larger  time  to  allow 
for  a  substantial  movement  of  the  boundary. 


Table  4.  Diffusion  Equation  (9a)  with  Cross-Derivatives 
and  Dirichlet  Boundary  Conditions  at  5=0, 

21 x  41  Grid 


T 

NSTEP 

NFE 

NJE 

0.05 

55 

77 

9 

0.10 

67 

90 

11 

0.15 

74 

98 

12 

0.20 

80 

105 

13 

0.25 

85 

110 

13 

0.  30 

89 

116 

13 

0.35 

93 

123 

!  14 

0.40 

96 

127 

14 

0.45 

99 

132 

14 

0.50 

103 

140 

14 

CPU  (sec) 


0.0026 

0.0055 

0.0076 

0.0102 

0.0126 

0.0126 

0.0150 

0.0150 

0.0150 

0.0150 


11.259 


16.911 
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For  our  final  case,  we  solve  equation  (9)  with  a 
nonuniformly  moving  boundary  and  the  correct  boundary 
condition  (9e) .  Recall  that  a  second  approximation  to 
the  Jacobian  is  used  with  this  boundary  condition.  The 
numerical  results  are  illustrated  in  Table  5.  The  number 
of  internal  iterations  is  now  as  much  as  65%  larger  than 
the  previous  case.  Also,  the  time  steps  are  smaller. 
Clearly,  the  boundary  approximation  to  the  Jacobian  is  a 
limiting  factor  on  computer  time.  In  the  future  we  will 
attempt  to  improve  this  approximation.  Nevertheless,  the 
computer  time  is  quite  reasonable  for  such  a  difficult 
problem. 


Table  5.  Numerical  Solution  of  Equation  (9), 
21 x  41  Grid 


T 

NSTEP 

NFE 

NJE 

Nil 

NQ 

At 

CPU  (sec) 

0.05 

70 

92 

19 

180 

4 

0.0022 

0.10 

86 

112 

■1 

232 

4 

0.0040 

0.15 

97 

131 

13 

288 

4 

0.0062 

0.20 

105 

143 

13 

323 

4 

0.0074 

0.25 

111 

157 

13 

365 

4 

0.0074 

0.30 

118 

178 

14 

428 

4 

0.0074 

0.35 

123 

189 

15 

463 

4 

0.0099 

0.40 

128 

204 

15 

517 

4 

0.0099 

0.45 

134 

220 

15 

574 

3 

0.0084 

0.50 

140 

231 

16 

612 

3 

0.0084 

27.97 

L 
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4.0  FUTURE  WORK 

To  date,  numerical  results  in  two  dimensions  have  been 
obtained  only  for  the  method  of  lines  algorithm,  and  nonuni¬ 
form  motion  of  the  oxide-silicon  boundary  has  been  handled 
numerically  by  the  translation-stretching  transformation. 

We  are  currently  programming  the  conformal  map  procedure 
outlined  in  Section  3.1.2  as  a  potentially  faster  way  to 
deal  with  moving,  nonrectangular  silicon  boundaries  typical 
of  the  formation  of  bird's  beak  structures  in  the  oxide  layer. 

In  the  near  term,  our  emphasis  will  be  on  evaluation  of 
the  method  of  lines  algorithm  against  realistic  test  cases, 
some  of  which  can  be  solved  to  known  accuracy  by  other 
techniques  (e.g.,  linear  diffusion  in  an  impenetrable  box). 

If  time  permits,  one  or  two  other  algorithms  which  are 
potentially  very  fast  will  be  compared  against  the  method 
of  lines. 

By  the  end  of  the  first  contract  year,  the  speed  and 
accuracy  of  the  basic  algorithm  should  be  well  characterized 
for  typical  process  steps.  The  relative  merits  of  different 
approaches  to  the  treatment  of  a  nonuniformly  moving  oxide- 
silicon  boundary  should  also  be  established. 

The  second  year  of  the  contract  will  be  devoted,  as 
planned,  to  seeking  algorithms  for  the  solution  of  more 
realistic,  physically  based  models  for  the  evolution  of 
VLSI  device  structures.  Our  objective  is  to  provide  the 
device  physicist  with  computational  analysis  tools  which 
can  with  tolerable  speed  and  accuracy  characterize  microscopic 
processes  involving  several  defect  species  migrating  under 
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the  influence  of  chemical,  mechanical,  and  electrical  driving 
forces  in  two  (and  possibly  three)  dimensions. 

This  work  will  also  provide  the  background  for  construct¬ 
ing  fast  algorithms  for  computer  modeling  of  complete  device 
electrical  characteristics  based  on  two-dimensional  dopant 
profiles.  While  the  coupled  system  of  carrier  transport  and 
Poisson's  equations  is  considerably  more  complex  than  the 
equations  describing  multispecies  diffusion,  there  is  a 
strong  analogy  with  the  equations  governing  steady-state 
fluid  flow  which  can  be  exploited. 
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